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1. Introduction 


Numerical inverse Laplace transform (NILT) methods are widely used in various scientific 
areas, especially for a solution of respective differential equations. In field of an electrical 
engineering many various approaches have been considered so far, but mostly for a single 
variable (1D NILT), see at least (Brancik, 1999, 2007b; Cohen, 2007; Valsa & Brancik, 1998; 
Wu at al., 2001) from plenty of papers. Much less attention was paid to multidimensional 
variable (nD NILT) methods, see e.g. (Hwang at al., 1983; Singhal at al., 1975), useful rather 
for more complicated electromagnetic systems. The 2D NILT methods, see e.g. (Brancik, 
2005, 2007a, 2007b; Hwang & Lu, 1999), can be applied for a transmission line analysis, or 
nD NILT methods, n 2 2, for a nonlinear circuits analysis, if relevant Laplace transforms are 
developed through a Volterra series expansion, see e.g. (Brancik, 2010a, 2010b, Karmakar, 
1980; Schetzen, 2006), to highlight at least a few applications. This paper is focused on the 
class of NILT methods based on complex Fourier series approximation, their error analysis, 
their effective algorithms development in a Matlab language, and after all, on their selected 
applications in field of electrical engineering to show practical usefulness of the algorithms. 


2. Multidimensional numerical inverse Laplace transform 


An n-dimensional Laplace transform of a real function f(t), with t = (t,...,4:) as a row vector 
of n real variables, is defined as (Hwang at al., 1983) 


oo 


bai n 
F(s)= f| -> f f()exp(-se")] Jdt, (1) 
ua fold 0 i-1 
where s = (51,...,51) and T means a transposition. Under an assumption | f(t) | < Mexp(at"), 
with M real positive and œ= (04,...,@,) being a minimal abscissa of convergence, and the nD 
Laplace transform F(s) defined on a region {s e C": Re[s] > aj, with c= (c1,...,c,) as an 
abscissa of convergence, and the inequality taken componentwise, the original function is 
given by an n-fold Bromwich integral 


1 Cytjoo ec, + joo 


Ga Í me Í Fsjexp(st Jas (2) 


CEJ Cage 
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In the papers (Bran¢ik, 2007a, 2007b, 2010b), it was shown for the 1D, 2D, and 3D cases, the 
rectangular method of a numerical integration leads to an approximate formula whose a 
relative error is adjustable, and corresponds to the complex Fourier series approximation of 
a respective dimension. The method has been generalized for an arbitrary dimension n in 
the recent work (Bran¢ik, 2011). 


2.1 Complex Fourier series approximation and limiting relative error 

Substituting s;=c;+ ja; into (2), and using a rectangular rule of the integration, namely 
@ = m Qi, and Q; = 27/7 as generalized frequency steps, with 7 forming a region of the 
solution te [0,7%) x...x [0,z,), an approximate formula is 


OPES Sma) 


M,=-99 m,=-% 


with s; = c; + jm;Q;, Vi. As is shown in (Brancik, 2011), a limiting relative error dy of (3) can 
be controlled by setting c = (c1,...,cn), defining paths of the integration in (2), namely 


Ci 


1 1 1, ô 
E e ta ye M, (4) 
"4 | ae] qon 
for i = 1,...,n, and while keeping the equalities n(c - @%) =...= (Cn - æn). The simplification 
in (4) is enabled due to small values dy considered in practice. The last equation is used for 
setting up parameters of the nD NILT method relating them to a limiting relative error dy 
required for practical computations. 


2.2 Practical computational methods 

It should be highlighted that the formula (4) is valid, and a relative error supposed is really 
achievable by the nD NILT (3), if infinite numbers of terms are used in the series. In practice, 
it cannot sure be fulfilled, but a suitable technique for accelerating a convergence of infinite 
series is usable, as is e.g. a quotient-difference (q-d) algorithm (Rutishauser, 1957). Besides, 
as has been already successfully used for cases of n < 3, the formula (3) can be rearranged to 
enable using FFT & IFFT algorithms for an effective computation. 


2.2.1 Partial ILTs evaluation technique 
The technique of practical evaluation of the n-fold infinite sum (3) follows the properties of 
the n-fold Bromwich integral (2), namely we can rearrange it into the form 


1 cy + joo 1 Cy + joo 1 C, +jo 
(ligt HS — --—— | F(s,,55,...,8, eds,» eds, leds , (5) 
AO) nil da eh OT a Sie 
or shortly 
f(titor rtn) z ia Ly ED [Eles] ` (6) 


Although the order of the integration may be arbitrary on principle, here the above one will 
be used for an explanation. Similarly, (3) can be rewritten as 
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Ti my =—% Ta m =—00 n m,=—% 


t eth z eo2'2 Z eon'n JM, Qant jm Qt jm Qyt 
f(tistar ta) = £ £ kes m a F(S1,897-++/Sy Je mening ag |17254212 | oJ IAL , (7) 


with s; = ci + jm;iQ;. If we define F, = F(s1,...,Sn-1,$n) and Fo = f(t1,...,tn-1,tn), then n consequential 
partial inversions are performed as 


Ly {F,} = Eal Sack 


12n-17 "n 


Lpa {En-a} = Ep-2(S1r ortn- , (8) 


Ly =F e A 


As is obvious we need to use a procedure able to make the inversion of Laplace transforms 
dependent on another n-1 parameters, complex in general. Let us denote arguments in (8) by 
Pi = (Pr---Pn-1,Pn). Then the ILT of the type 


C; +jo 


Falp) = GHE) m | Eppe" ds; (9) 


c-j 


can be used n times, i = n,n-1...,1, to evaluate (8), with pn = (S1,...,Sn-1/Sn), Pred = (St--/Sn-trtn) yo 
pi = (S1--tn-ttn), and po = (t1..,tnatn), while pj =s; for j <i, and p;=t; otherwise. A further 
technique is based on demand to find the solution on a whole region of discrete points. 
Then, taking into account tix = kT; in (9), with T; as the sampling periods in the original 
domain, we can write an approximate formula 


c;kT; œ 


Daye ee (10) 


1 m=—°eo 


e 


PaPa = 


i = n,n-1..,1, and with Q;= 27/7 substituted. As follows from the error analysis (Bran¢ik, 
2011) a relative error is predictable on the region Oer = [0,%) *...* [0,7%). For k = 0,1,...,.Mi-1, 
i = 1,...,n, a maximum reachable region is Omax = [0,(Mi-1)T1] *...x [0,(M,-1)T,,]. Thus, to meet 
the necessary condition Omax © Oer, we can set up fittingly 7 = M;iT, i =1,....n. In practice, a 
region of the calculation is chosen to be Oea = [0,teai] *...% [O,tncail, with tica = (Mi/2-1)Ti, 
i = 1,....n, to provide certain margins. 


2.2.2 FFT, IFFT, and quotient-difference algorithms utilization 

As is shown in (Brancik, 2007a, 2010c), the discretized formula (10) can be evaluated by the 
FFT and IFFT algorithms, in conjunction with the quotient-difference (q-d) algorithm for 
accelerating convergence of the residual infinite series, see following procedures. 

To explain it in more detail, let us consider an r-th cycle in gaining the original function via 
(9), i.e. F,_,(p,4)=L;,'{E,(p,)}. For its discretized version (10) we have 


E (s kT t ) = eorkT, ŞS Ë (s Pee 20 i j2amkT, Jr, 11 
a kt E a i preeeren m r\Opreresly JM— peeves aye x ( ) 
r è m=—% r 


The above stated formula can be decomposed and expressed also as 
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ec kT, M,- co 
Ë lsi kl euch) = S BCom 4 yess S3 m > Gm — FO |, (12) 
T, m=0 m=0 n=0 m=0 
where individual terms are defined as 
M, = 26, K, integer , 
FE = È (5,,...,¢, + JMT |T, p.t); (13) 
=(£m) — FEM, 4» +m) 
r r , 


when ge = etl2ak =] , Vk , has been considered, and 7 = M,T,. 

As is evident the first and the third finite sum of (12) can be evaluated via the FFT and IFFT 
algorithms, respectively, while 2P+1 terms from the infinite sums are used as the input data 
in the quotient-difference algorithm (Macdonald, 1964; McCabe, 1983; Rutishauser, 1957). 
We can replace the above infinite power series by a continued fraction as 


DG ath, = do/(1+ dzax /(1+ + + dopZay)) Wk, (14) 


which gives much more accurate result than the original sum truncated on 2P+1 terms only. 
The q-d algorithm process can be explained based on a lozenge diagram shown in Fig. 1. 


20 
0 
qi” 
el) el) 
qf” gP 
0) at 2) 
qi” qs 
e) el?) 
3 
4") 
e 


Fig. 1. Quotient-difference algorithm lozenge diagram 


The first two columns are formed as 


e® =0, i=0,...,2P, 
1 ~ti (15) 
gi) = GROG, i=0,...,2P-1, 
while the successive columns are given by the rules 
(i) = gl Wa gt) j= _ 
es =q" +e, 1=0,-+,2P—2)], for j=1,..,P, 
j qj” Jj teja (16) 


gP = qt! o i=0,---,2P-2j-1, for j=2,...,P. 
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Then, the coefficients dm, m = 0,...,2P, in (14) are given by 


=F a, = 0 = en 
dy= GO, dja, dj-e, j=1..P. (17) 
For practical computations, however, the recursive formulae stated below are more effective 
to be used (DeHoog et al., 1982). They are of the forms 


An(Zsx) = Am1(Z4x) F An Z+¢Am—2 (ZK) , Bi (Z44) = Bw (Zsje) +t Ain Z-¢Bm_2 (Zak) , (18) 


for m =1,...,2P, Vk, with the initial values A. = 0, B = 1, Ao = do, and Bo = 1. Then, instead of 
the continued fraction (14), we can write 


Ps Gh = Aple Bolir) VE (19) 
m=0 


The q-d algorithm is a very efficient tool just for a power series convergence acceleration, 
here enabling (7) to achieve a relative error near its theoretical value defined by (4), see the 
following examples. 


2.3 Matlab listings and experimental errors evaluation 

In this part experimental verifications of the nD NILT theory above will first be presented, 
for one to three dimensional cases, i.e. n < 3. For such dimensions the Matlab functions have 
been developed and errors stated on a basis of some sample images with known originals. 
The Matlab listings of basic versions of the NILT functions are provided, together with 
examples of their right calling. Another Matlab listings will be discussed in more detail later, 
in the chapter with practical applications. 


2.3.1 One-dimensional NILT 
In case of the 1D inverse LT, a well-known Bromwich integral results from (2), namely 


1 c+ joo 
fH=s | Feta, (20) 
20 j ire 
where indexes 1 were omitted. By using the theory above a path of the numerical 
integration is stated according to (4), leading to 


eee@=inlit=—— \=e4 lies lee- he ,. (21) 
T 1+ dy T Ôm T 


In contrast to most other approaches, the 1D NILT method described here enables to treat 
complex images resulting in complex originals as no real or imaginary parts are extracted 


during an evaluation process. It can be useful in some special applications, not only in the 
electrical engineering. We can śhow it on a simple transform pair 


F(s)= = tists B f= e= coso jsinat. (22) 
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Of course, when preprocessing the transform to arrange it to a Cartesian form, as is shown 
on the right sides in (22), the result could by get by inverting the real and imaginary parts 
separately, by using an arbitrary NILT method. Here, however, no symbolic manipulations 
are needed in advance, and F(s) enters the NILT function in its basic form as a whole. 

A Matlab language listing is shown in Tab. 1, where the relative error needed is marked by 
Er and is subject to a change if necessary, similarly as the minimal abscissa of convergence 
(exponential order) œ, alfa, numbers of points for the resultant solution, M, and for the q-d 
algorithm, P. If only real transforms F(s) are considered the bottom line in the listing can be 
inactivated. The NILT function is called from a command line as follows: niltc('F',tm) ; 
where F is a name of another function in which the F(s) is defined, and tm marks an upper 
limit of the original variable t. In our case, and for œ= 27, this function can have a form 


function f=expc(s) 
f=1./(s-2*pi*j); 


———=— 1D NILT for complex arguments - basic version 
TA AE based on FFT/IFFT/q-d, by L. Brančík 
function [ft,t]=niltc(F,tm); 
alfa=0; M=256; P=3; Er=le-10; 
N=2*M; qd=2*P+1; 
t=linspace(0,tm,M); NT=2*tm*N/ (N-2); omega=2*pi/NT; 
c=alfatlog(1+1/Er)/NT; s=c-i*omega* (0:N+qd-1) ; 
Fs(1,:)=feval(F,s); Fs(2,:)=feval(F,conj(s)); 
ft(1,:)=fft (Fs(1,1:N)); £t(2,:)=N*ifft (Fs(2,1:N)); 
ft=ft(:,1:M); D=zeros(2,qd); E=D; 
Q=Fs(:,N+2:N+qd) ./Fs(:,N+1:N+qd-1) ; 
D(:,1)=Fs(:,N+1); D(:,2)=-Q(:,1); 
for r=2:2:qd-1 
w=qd-r; 
E(:,l:w)=Q0(:,2:wtl)-Q(:,1l:w)t+E(:,2:wtl); D(:,rt+1)=-E(:,1); 
if r>2 
Q(:,1l:w-1)=0(:,2:w).*E(:,2:w)./E(:,1l:w-1); 
D(:,r)=-Q(:,1); 
end 
end 
A2=zeros(2,M); B2=ones(2,M); Al=repmat (D(:,1),[1,M]); B1=B2; 
zl=exp (-i*omega*t); z=[z1;conj(z1)]; 
for n=2:qd 
Dn=repmat (D(:,n),[1,M]); 
A=A1+Dn.*z.*A2; B=B1+Dn.*z.*B2; A2=Al; B2=Bl; A1l=A; B1=B; 
end 
ft=fttA./B; ft=sum(ft)-Fs(2,1); ft=exp(c*t)/NT.*ft; 
ft (1) =2*ft (1); 
figure; plot (t,real(ft) 
figure; plot(t,imag(ft)); % optional 


Table 1. Matlab listing of 1D NILT method accepting complex arguments 


As is obvious, the Laplace transform must be defined to enable Matlab array processing, i.e. 
element-by-element array operators have to be used. Thus, the calling our function can look 
like niltc('expc', 4); if the function is saved under the same name, expc, or it is placed 
inside the M-file with own NILT function (Tab. 1), following always its body. Alternatively, 
the calling can look like [ft,t]=niltc('F',tm)j; if respective variables in the brackets 
are to be saved in the memory after the function finishes. 
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Graphical results and corresponding errors are shown in Fig. 2. Because the originals are 
bounded by values +1, and w= 0, we can see the errors satisfy (21) very well (dy = 10-10 was 
considered), excluding only beginning of the interval. 


Real part of the original Imaginary part of the original 
1 T T 


Fig. 2. Numerical inversion leading to complex original f(t) = exp(jat) 


Another test functions are considered in Tab. 2, with numerical results shown in Fig. 3. As is 
again obvious from Fig. 3 the relative errors satisfy theoretical expectations, with an 
exception of vicinities of discontinuities. 


Table 2. Test Laplace transforms for errors evaluation 


2.3.2 Two-dimensional NILT 
In case of the 2D inverse LT, a two-fold Bromwich integral results from (2), namely 


1 Ci + jo C, + joo 
Ie)? = f f F(s,,8,)e"'*? ds,dsy , (23) 


c= Joe Cy — Joe 


and by using the theory above the paths of numerical integrations are stated based on (4) as 


1 1 1, ô 
¢; = ar, -—In| 1- ——— |=a,-—In—,, i=1,2. 24 
a | a] a2 e 
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ORIGINALS r ERRORS ORIGINALS i ERRORS 
r r F T T r 1 r r 0 r r 


Fig. 3. Computed originals and errors for test Laplace transforms in Tab. 2 


A Matlab language listing is shown in Tab. 3, with all the parameters denoted by similar 
way as in the previous case. 

Nevertheless, the Laplace transform variables and the original variables have changed in 
this listing as sı > p, s2 > q, and tı — x, t2 > y, respectively, which simplified a writing. The 
parameters are then indexed in compliance with these new notations. Besides, numbers of 
points used to plot three-dimensional graphical results are set by xp1 and yp1. 

For the same reasons as at the 1D NILT, the 2D NILT method discussed here enables to treat 
complex images of two variables resulting in complex originals. We will Show it on a simple 
transform pair 


F(s,,52) = > f (tt) = gate ata) . (25) 


(sı - jæ )(s2 — ja) 
After rearranging the above equation, we can also write 
$85 — l s, + @s 
F(s,,8)) = 2 = an zT] A I 7n” 
(s; + af )(s3 +a) (si +a \(s3 +@5) f (26) 
f(t tr) = cos(@t, + @t,) + jsin( aft, +t) 


The 2D NILT function is called from a command line as follows: nilt2c('F',xm, ym); 
where F is a name of another function in which the F(p,q) is defined, and xm and ym mark 
upper limits of the original variables x and y. In our case, and for @ = @ = 22, this function 
can have a form 


function f=exp2c (p,q) 
f=1./(p-2*pir*j) ./(q-2*pi*j); 


and its calling can look like nilt2c('exp2c', 3,3); with graphical results in Fig. 4. As 
the originals are bounded by values +1, and a@ = œ = 0, we can see the errors satisfy (21) 
very well (dv = 10° was considered), excluding beginnings of the 2D region. 
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zess 2D NILT based on partial inversions, by L. Brančík 
function fxy=nilt2c(F,xm, ym) ; 
alfax=0; alfay=0; Mx=256; My=256; P=3; Er=le-8; % adjustable 
xpl=64; ypl=64; % adjustable 
Nx=2*Mx; Ny=2*My; qd=2*P+1; Ke=log(1-1/sqrt (1+Er) ); 
nx=2*xm*Nx/ (Nx-2); ny=2*ym*Ny/ (Ny-2) ; 
omegax=2*pi/nx; omegay=2*pi/ny; sigx=alfax—Ke/nx; 
sigy=alfay-Ke/ny; qdl=qd-1; Nxw=Nx+qd1; Nyw=Ny+qd1; 
Asigx=sigx-i*omegax* (0:Nxw); Asigy=sigy-—i*omegay* (0:Nyw) ; 
Asigx2=cat (2,Asigx, conj (Asigx) ); 
rx=[1:Mx/xpl:Mx, Mx]; ry=[1:My/ypl:My,My]; 
x=linspace (0, xm, Mx); y=linspace(0,ym,My); x=x(rx); y=y(ry); 
[q, p]=meshgrid (Asigy,Asigx2); Fpq(:,:,1)=feval(F,p,q); 
[q, p]=meshgrid (conj (Asigy),Asigx2); Fpq(:,:,2)=feval(F,p,q); 
Fpyp=Pnilt (Fpq, Ny, ry, qd, Yy, ny, omegay, sigy); % Pnilt to get F(p,y) 
Fpy (:, :,1)=Fpyp(1:Nxw+1,:).'; 
Fpy(:,:,2)=Fpyp (Nxw+2:2*Nxw+2,:).'}; 
fxy=Pnilt (Fpy,Nx, rx, qd, x, nx, omegax, Sigx) ; 5 lt to get f(x,y) 
figure; mesh (x,y, real (fxy)); 
figure; mesh (x,y, imag (fxy)); % optional 
=== PARTIAL NILT based on FFT/IFFT/Q-D, by L.Branéik 
function fx=Pnilt (Fq,N, grid, qd, xy, nxy, omega, C); 
fx(:,:,1)=fft(Fq(:,:,1),N,2); £x(:,:,2)=N*ifft (Fq(:,:,2),N,2); 
fx=fx(:,grid,:); delv=size(Fq,1); delxy=length (xy); 
d=zeros(delv,qd,2); e=d; g=Fq(:,N+2:N+qd,:)./Fq(:,N+1:Nt+qd-1,: 
d(:,1,:)=Fq(:,Nt+1,:); d(:,2,:)=-q(:,1,:); 
for r=2:2:qd-1 
w=qd-r; e(:,l:iw,:)=q(:,2:wtl,:)-q(:,liw,:)+e(:,2:wtl,:); 
ad(:,xrt1, 2:)=-e(4,1; 2); 
if eZ 
q(:,l:w-1,:)=q(:,2:w,:).*e(:,2:w,:)./e(:,l:iw-1,:); 
a(:,r,:)=-q(:,1,:);7 
end 
end 
A2=zeros (delv,delxy,2); B2=ones (delv,delxy, 2); 
Al=repmat (d(:,1,:),[1,delxy,1]); B1=B2; 
z1(1,:,1)=exp(-i*omega*xy); z1(1,:,2)=conj(z1(1,:,1)); 
z=repmat (z1, [delv,1]); 
for n=2:qd 
Dn=repmat (d(:,n,:),[1,delxy,1]); 
A=A1+Dn.*z.*A2; B=B1+Dn.*z.*B2; A2=Al; B2=Bl; A1l=A; B1=B; 
end 
fx=fxtA./B; fx=sum(fx,3)-repmat (Fq(:,1),[1,delxy]); 
fx=repmat (exp (c*xy) /nxy, [delv,1]).*fx; fx(:,1)=2*fx(:,1); 


Table 3. Matlab listing of 2D NILT based on partial inversions 


Another simple example shows a shifted 2D unit step, with different shifts along the axis. A 
corresponding transform pair is 


exp(—2s, —s 
F (5,8) = EPES) 5 (ty ts) =U -2-1), (27) 
1°2 


In this case, a displaying imaginary part gives a zero plane, and the respective line in the 2D 
NILT function can be inactivated. The graphical results are depicted in Fig. 5, including an 
absolute error. The respective Matlab function can be of a form 
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function f=step2 (p,q) 
f=exp (-2*p-q) ./p./q; 


and called as nilt2c('step2', 4,4) ;, with the results theoretically expected. 


Real part of the original Imaginary part of the original 


Imif(t,.t,)] 


Fig. 4. Numerical inversion leading to complex original f(t1,t2) = exp(j@(ti+t)) 


Shifted 2D unit step Absolute error for the 2D unit step 


i Hii 2 
it Whi] iili O 
Hh FN iH : = 
| a | i i i (Hil Hi Ww 

intl 


Fig. 5. Numerical inversion leading to shifted 2D unit step f(t,t2) = 1(t—-2,t2-1) 


2.3.3 Three-dimensional NILT 
In case of the 3D inverse LT, a three-fold Bromwich integral results from (2), namely 


. + joo Cy + joo C3 + joo 
Taa A [J Elsisas5)e 22794 ds ds ds; (28) 


— Jor Cy — jæ C3 — joo 


and by using the theory above the paths of numerical integrations are stated based on (4) as 
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ae a dede a 

Ci = nf: E 7 In 3 i=1,2,3. (29) 
Here only experimental results will be shown to verify an accuracy of the method. A Matlab 
language listing looks similarly like for the 2D NILT case, but the partial NILT subfunction 
is called once more, and respective arrays dimensions are enlarged. Original functions 
corresponding to 3D Laplace transforms cannot be displayed graphically as a whole, of 
course. However, for one variable chosen as constant, it is posssible to display three 
respective two-dimensional cuts. It will be demonstrated on the example of 3D shifted unit 
step, with a Laplace transform pair 


exp(—s, — 2s, —3s3) 


e (t-l, -—2,t,-3), (30) 
515953 


with different values of shifts along respective coordinates so that correctness of results can 
easily be identified, see Fig. 6. Errors again correspond to theoretically expected ones. 


Unit-step cut I(t, ,t,), t, = const. Unit-step cut 1(t,,t,), t, = const. Unit-step cut I(tyt,), t, = const. 


3 


Fig. 6. Numerical inversion leading to shifted 3D unit step f(t,t2,t3) = 1(t-1,t2-2,t3-3) 


3. Application of NILT algorithms to electrical engineering simulation 


In this chapter some examples of the application of the NILT algorithms developed relating 
to problems of electrical engineering simulation are presented. First, the 1D NILT method is 
applied for the solution of transient phenomena in linear electrical circuits with both 
lumped and distributed parameters. This well-known approach is usable wherever linear 
ordinary differential equations (ODE) are transformed into algebraic ones so that an inverse 
Laplace transform can be considered. Then the 2D NILT method is utilized to solve transient 
phenomena on transmission lines (TL) after relevant telegraphic equations (a type of partial 
differential equations (PDE)) are transformed into algebraic ones by a 2D Laplace transform. 
In this way voltage and/or current distributions along the TL wires can be determined in a 
single calculation step. Finally, the utilization of the 1D to 3D NILTs to weakly nonlinear 
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electrical circuits solution is discussed. In this case the relevant nonlinear ODEs describing 
the circuit are expanded into Volterra series which respective NILTs are applied on. 


3.1 One-dimensional NILT algorithm application 

3.1.1 Preliminary example based on lumped parameter circuit 

A simple example demonstrating the application of the basic 1D NILT algorithm in Tab. 1 is 
shown in Fig. 7. This really initiatory linear electrical circuit was chosen with an intention to 
be also considered later, in chapter 3.3.1, as a nonlinear circuit, with G2 being a nonlinear 
element. In this way one will be able to compare results and make some conclusions. 


Fig. 7. Linear reactive electrical circuit of the 1st order 


Denoting G = G1 + G, the 1s-order linear ODE has a form 


cl) 
dt 


+Gv(t)=i(t), (31) 
with a Laplace-domain solution 


V(s)= I(s) + Cv(0) (32) 
Gt+sC 
with an initial condition v(0). Even if the above circuit is very simple a finding time-domain 
solution could be rather work-intensive if the circuit is excited from some non-trivial input 
current waveform. A few basic examples are given in Tab. 4, specially the first one results in 
a transient characteristic of the circuit. 


3 4 


I)sin(2zt)1(t) | Ip cos(2zt)1(t) 


2719 slo 
3? +47 s? +47 


Table 4. Exciting current source waveforms and their Laplace transforms 


For the above examples, of course, time-domain analytical solutions can be found e.g. based 
on a Heaviside formula. The 1D NILT function graphical results, under a condition v(0) = 0, 
and considering values C = 1mF, G1 = G2 = 10mS, and Ip = 1mA, are shown in Fig. 8. 
The above waveforms can be got by either successive application of a basic version of the 1D 
NILT method according to Tab. 1, or a generalized 1D NILT function, its vector version, can 
be used to process all the computations in parallel. This function is shown in Tab. 5. 
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x10° Exciting current waveforms Voltage responses for exciting currents 
1.5 T l l j 0.05 


tis l tis] 


=SS=s= 1D NILT for complex arguments - vector version ------- 
Seend based on FFT/IFFT/q-d, by L. Brančík ------------ 

function [ft,t]=niltcv (F,tm,depict); 
alfa=0; M=256; P=3; Er=le-10; adjustable 
N=2*M; qd=2*P+1; t=linspace(0,tm,M); NT=2*tm*N/ (N-2); 
omega=2*pi/NT; 
c=alfatlog(1+1/Er)/NT; s=c-i*omega* (0:N+qd-1) ; 
Fs(:,:,1)=feval(F,s); Fs(:,:,2)=feval(F,conj(s)); lv=size(Fs,1); 
ft (z; L)y=ffFt (Fs (i; 271), N2); fE Ct, 2, 2)SN* 2 £EC (FS (2; 27-2) 7 Np 2) 7 
ft=ft(:,1:M,:); 
D=zeros(lv,qd,2); E=D; Q=Fs(:,N+2:N+qd,:)./Fs(:,N+1:N+qd-1,:); 
D(:,1,:)=Fs(:,Nt1,:)% D(:,2,:)=-Q(:,1,:); 
for r=2:2:qd-1 

w=qd-r; 

E(:,1l:w,:)= 2)-O (2, 1liwy,:)4+E(: 

D(a EFL) Bp EA 

LE 2 

Q(:,1:w-1,:)=Q(: 

D(:,r,:)=-Q(:,1,: 

end 
end 
A2=zeros(lv,M,2); B2=ones(lv,M,2); Al=repmat (D(:,1,:),[1,M,1]); 
B1=B2; zl=repmat (exp (-i*omega*t),[lv,1]); z=cat(3,z21,conj(z1)); 
for n=2:qd 

Dn=repmat (D(:,n,:),[1,M,1]); 

A=A1+Dn.*z.*A2; B=B1+Dn.*z.*B2; A2=Al; B2=Bl; A1=A; B1=B; 
end 
ft=ftt+A./B; ft=sum(ft,3)-repmat (Fs(:,1 
ft=repmat (exp(c*t)/NT, [lv,1]) .*ft; ft(:,1) 
switch depict 

case 'pl', plottli(t,ft); case 'p2', plott2(t,ft); 

case 'p3', plott3(t,ft); otherwise display('Invalid Plot'); 


72),[1,M,1]); 
=2*ft(:,1); 


end 


Table 5. Matlab listing of vector version of 1D NILT method 


Here one more parameter depict is used to define a method of plotting individual items 
from a set of originals. The 1D NILT function is called as niltcv('F',tm, 'depict"); 
where 'depict' isa text string 'p1', 'p2"', or 'p3"', see Tab. 6 for more details. 
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--- Plotting functions called by 1D NILT, vector version ---- 
——S====—S-= Multiple plotting into single figur 
function plotti1(t, ft) 
figure; plot(t,real(ft)); grid on; 
figure; plot(t,imag(ft)); grid on; optional 
% Plotting into separate figures 
function plott2(t,ft) 
for k=1:size(ft,1) 
figure; plot (t,real (ft ( 
figure; plot (t,imag(ft ( 


)) 
)) 


grid on; 
grid on; % optional 


k,:) 
Ket) 


r 
r 


Plotting into 3D graphs 

function plott3(t, ft) 

global x; must be global in F 
m=length(t); tgr=[1:m/64:m,m]; time points chosen 
figure; mesh(t(tgr),x,real(ft(: 
figure; mesh(t(tgr),x,imag(ft(: ; optional 


Table 6. Matlab listing of plotting functions for vector version of 1D NILT method 


To get e.g. the right part of Fig. 8, that is the voltage responses of the circuit in Fig. 7, the 
calling the 1D NILT function looks like niltcv('V4',1,'p1'); where V4 denotes a 
name of the function defining individual responses as follows: 


function f=V4(s) 
TO=le-3; C=le-3; G=2e-2; 
£(1,:)=10./s./ (G+s*C); 
£(2,:)=10./(s+5)./(G+s*C); 
2) 
:) 


f (3, :)=2*pi*I0./(s.^2+4*pi^2)./(G+s*C); 
f (4,:)=s.*I0./(s.^2+4*pi^2)./(G+s*C); 


In this case the lines causing the imaginary parts plotting can be inactivated. The remaining 
plotting functions will be explained in the next chapter. 


3.1.2 Application for transmission line simulation 
Here, the 1D NILT algorithms will be used to simulate voltage and/or current distributions 
along transmission lines (TL), as shown on a Laplace-domain TL model in Fig. 9. As is well 
known, this model results from the application of a Laplace transform, with respect to time, 
on a pair of partial differential equations (telegraphic) of the form 

dv(t,x) p. dift, x) di(t, x) 

T Roi(t,x) + Lo at |” = 

with Ro, Lo, Go, and Co as per-unit-length (p.-u.-l.) parameters, being constant for uniform 
TLs, and with a length 1. 
When considering zero initial voltage and current distributions, v(0,x) = 0 and i(0,x) = 0, and 
incorporating boundary conditions, we get the Laplace-domain solution in the forms 


= Gyo(t,x) +c, EA) (33) 


-y(s)x —/(s)[21-x] 
Vsa) = Vils) E 6 ee (34) 


i(s)+Z.(8)_ 1-py(s)p2(s)e77™ 


I(s,x) = V;(s) l Eao ne 
Zi(8)+ Z8) 1- p (S)pa(9 
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hies I(s,x) h(s} 


ZA) 


Fig. 9. Laplace-domain model of transmission line with linear terminations 


where Z,(s) and %s) are a characteristic impedance and a propagation constant, respectively, 


R,+sL 
Z,(s)= ==  , = (Ra +sLo (Ga +sCo) , 36 
c(S) G,+sC, x(s) ( ofS o)( ots 0) (36) 


and p(s) and p(s) are reflection coefficients at the TL beginning and end, respectively, 


Zi(s) —Z-(s) 
Z;(s) + Z,(s) 


- 22(8)-Z-(s) 
>  Pr(s)= Zale). (37) 


p(s) = 
In a general case of lossy TLs, the time-domain solutions cannot be found by an analytical 
method, thus the only way is to use some numerical technique. 
As an example, let us consider the TL of a length / = 1m, with p.-u.-l. parameters Ro = 1mQ, 
Lo = 600nH, Go = 2mS, and Co = 80pF, terminated by resistive elements Z; = 10Q, Z2 = 1kQ, 
and excited by the voltage source waveform v;(t) = sin2(at/2 10°), 0 < t <2 10°, and v;(t) = 0, 
otherwise, with the Laplace transform 


E 2m’ [1 —exp(-2- 10s) | 


p s| (2:10°s)? +47” | os 


The Fig. 10 shows time dependences at the beginning, the centre, and the end of the TL, 
while the 1D NILT is called as niltcv('Vs',4e-8,'pl1'); where the function Vs is 
defined as 


function f=Vs(s) 

Ro=le-3; Lo=600e-9; Go=2e-3; Co=80e-12; 

Zi=10; Z2=1e3; 

Vi=2*pi*2* (1-exp (-2e-9*s))./s./((2e-9*8) .*2+4*pir2) ; 

Z=Ro+s*Lo; Y=Gots*Co; Zc=sqrt(Z./Y); gam=sqrt(Z.*Y); 

rol=(Zi-Zc) ./(ZitZc); ro2=(Z2-Zc) ./(22+2Zc); 

Ks=Vi./(ZitZc) ./(1-rol.*ro2.*exp (-2*gam*1) ); 

for k=1:length (x) 
f(k,:)=Ks.*Zc.* (exp (-gam*x (k) )+ro2.*exp (-gam* (2*1-x(k)))); 


end 


Similarly, current waveforms can be computed by the above function slightly modified 
according to (35). Both waveforms are depicted in Fig. 10. 

Finally, it will be shown, how to obtain three-dimensional graphical results representing 
voltage and current distributions along the TL. Besides a possibitity to use the for cycle, as 
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shown in the function Vs above, another method based on 3D arrays will be applied, see the 
function Vsx below: 


Voltage waveforms on TL Current waveforms on TL 


= 


T 
aa ee 
pane 


Fig. 10. Numerically computed TL voltage and current waveforms 


function f=Vsx (s) 

global x; 

1=1; 

Ro=le-3; Lo=600e-9; Go=2e-3; Co=80e-12; 

Zi=10; Z2=1e3; 

x=linspace(0,1,65); % 65 points along TL chosen 
[S,X]=meshgrid(s,x); 

Vi=2*pi*2* (1-exp (-2e-9*S)) ./S./((2e-9*S) .*2+4*pir2) ; 
Z=Ro+S*Lo; Y=GotS*Co; Zc=sqrt(Z./Y); gam=sqrt(Z.*Y); 
rol=(Zi-Zc) ./(ZitZc); ro2=(Z2-Zc) ./(2Z2+2Zc) ; 
Ks=Vi./(ZitZc) ./(1-rol.*ro2.*exp (-2*gam*1) ); 
f=Ks.*Zc.* (exp (-gam.*X)+ro2.*exp (-gam.* (2*1-X))); 


In this case, the 1D NILT algorithm in Tab. 5 is called as niltcv('Vsx',2e-8, 'p3'); that 
is the plott3 function is used for the plotting, see Tab. 6, and a time limit is half of that in Fig. 
10 to get well-observable results. Again, the current distributions can be gained via the above 
function slightly modified according to (35). Both 3D graphs are depicted in Fig. 11. 


Voltage distribution along the TL wire Current distribution along the TL wire 


0.015 


v(t,x) [V] 


t[s] x [m] 


Fig. 11. Numerically computed TL voltage and current distributions 


3.2 Two-dimensional NILT algorithm application 
A two-dimensional Laplace transform can generally be used for the solution of linear partial 


differencial equations with two variables. The advantage is that we get completely algebraic 
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equations leading to much easier solution in the Laplace domain. A final step in the solution 
is then the utilization of the 2D NILT algorithm to get results in the original domain. Such a 
possibility will be shown on the example of telegraphic equations describing transmission 
lines, and results will be compared with the 1D NILT approach. 


3.2.1 Application for transmission line simulation 

Herein, rather less conventional approach for the simulation of voltage and/or current 
distributions along the TLs will be discussed. As is obvious from the telegraphic equations 
(33) they can be transformed not only with respect to the time t, which was matter of the 
previous paragraph, but also with respect to the geometrical coordinate x to get completely 
algebraic equations. After performing such the Laplace transforms, incorporating boundary 
conditions given by the terminating circuits, and considering again zero initial voltage and 
current distributions, v(0,x) = 0 and i(0,x) = 0, we get (Valsa & Bran¢ik, 1998b) 


_ 1Vi(s)— Xs)Zc(s)hi(s) (39) 


1(s,q)= 9). _, (40) 


where V;(s) = V(s,0) and I;(s) = I(s,0) are given by (34) and (35), respectively, see also Fig. 9. 
Thus the 2D NILT function according to Tab. 3 can be called as nilt2c('Vsq',2e-8,1); 
leading to the same graphical results as are shown in Fig. 11. The function Vsq can be of the 
form as stated below. The current distribution is obtained via the same function slightly 
modified according to (40). 


function f=Vsq(s,q) 

1=1; Zi=10; 2Z2=1e3; 

Ro=le-3; Lo=600e-9; Go=2e-3; Co=80e-12; 

Vi=2*pi*2* (1-exp (-2e-9*s))./s./((2e-9*8) .*2+4*pir2) ; 
Z=Ro+s*Lo; Y=Gots*Co; Zc=sqrt(Z./Y); gam=sqrt(Z.*Y); 
rol=(Zi-Zc) ./(ZitZc); ro2=(Z2-Zc) ./(22+2Zc); 
Ks=Vi./(ZitZc) ./(1-rol.*ro2.*exp (-2*gam*1) ); 
V1l=Ks.*Zc.%* (1+ro2.*exp (-2*gam*1)); 

I1=Ks.* (1-ro2.*exp(-2*gam*1) ); 

f=(q.*V1-Zc.*gam.*I1) ./(q.*2-gam.%2) ; 


One can notice an interesting thing, namely getting both voltage and current graphs by a 
single computation step. It is enabled by putting together the voltage and current transforms 
forming respectively real and imaginary parts of an artificial complex transform, and letting 
active the program command for the plotting the imaginary part of the original function, see 
Tab. 3. In our example, if the bottom line in the Vsq function is changed to 


f=((q.*V1-Zc.*gam. *1I1)+4* (q.*I1l-gam./Zc.*V1))./(q.*2-gam.%*2) ; 
then both graphs in Fig. 11 are obtained simultaneously. The same possibility exists for the 


1D NILT functions discussed earlier. There is no obvious physical meaning of such articifial 
complex transfoms, it is only a formal tool for inverting two transforms in parallel instead. 
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3.3 Multidimensional LT in nonlinear electrical circuits simulation 

As is known some classes of nonlinear systems can be described through a Volterra series 
expansion, accurately enough from the practical point of view, when a response v(t) to a 
stimulus i(t) has a form (Schetzen, 2006) 


H(t) = Lott) (41) 


where the terms of the infinite sum are 


a n 


D,(t)= f pa | a(t t-te) I EE (42) 


—co —oo k=1 


with hn(%,T2,...,%) as an n-th order Volterra kernel, called as a nonlinear impulse response 
as well. The Fig. 12 shows these equations in their graphical form. 

By introducing new variables, t = h =...= t, = t, and by using the n-dimensional Laplace 
transform (1), the n-fold convolution integral (42) leads to a Laplace domain response 


Va (S1827 tSn) = Hp(SirS27 Sa) [ES 7 (43) 
k=1 


with H,,(s1,S2,...,$n) aS a nonlinear transfer function. 


Fig. 12. Nonlinear system response via Volterra series expansion 


A few methods are at disposal to find the transfer function for a given nonlinear system, like 
a harmonic input method, e.g. (Bussgang at al., 1974; Karmakar, 1980). Further procedure is 
usually based on the association of variables, (J. Debnath & N.C. Debnath, 1991; Reddy & 
Jagan, 1971), transforming V,,(s1,Sz,...,5) into the function of a single variable V,,(s), and 
enabling to use a one-dimensional ILT to get the required terms v;(t) in (41). In contrast to 
this procedure, it is also possible to determine the terms 0;(t1,t2,...,tn) by the use of the n- 
dimensional ILT, considering h = t = ...= t, = t in the result as a final step. That is why the 
above NILT procedures can be adapted in this respect being able to serve as a tool for the 
nonlinear circuits transient simulation. 
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3.3.1 Utilization of 1D to 3D NILTs for nearly nonlinear circuits 

The utilization of the NILT methods developed, up to three-dimensional case, will be shown 
on the solution of a nearly nonlinear circuit in Fig. 13. As can be observed this is just Fig. 7 
modified to introduce a nonlinearity via G7 conductance. 


Fig. 13. Electrical circuit with nonlinear resistive element G2 


Assuming a square nonlinearity, a circuit equation is 


c dv(t) 
dt 


+G,0(t) + Gov" (t) = i(t). (44) 


By using the harmonic input method, and limiting the solution on the first three terms only, 
we get the nonlinear transfer functions for (43) as 


Hi(s)= (sC +G) (45) 
H3 (81,82) = —GyH,(s,) Hy ($2 )Ay (sı +82) (46) 


G 
13 (5/82/83) = -= THiS: Ha (62-53) + H,(s)H5(s,,83) + Hy (s3)Hp(s,-52)] Hy (s; +82 +53) - (47) 
Let us use an exciting current and its Laplace transform as 


i(t)=Ije l(t) = I(s)= E (48) 
s+a 
a20. The substitution (45) - (48) into (43) gives us respective Laplace-domain responses 
which will undergo the 1D, 2D and 3D NILT algorithms, respectively. We can write 


v(t) = v(t) + v(t) + 03(£) = v(t), 4 + valita) e-t + vs(tirtarts)l op apt T 


| i i , (49) 
=L;'[V, ED, +L [vers] a +Ly [V3(s1,52/83)]| 


t=t,=t,=¢ 


with L;'[.] as a k-dimensional ILT. Thereby, a time-consuming association of variables can 


be omitted, e.g. (Wambacq & Sansen, 2010). Individual terms v(t) are depicted in Fig. 14, for 
values agreeing with the linear circuit version in Fig. 7. The current i(t) is defined by a = 0 (a 
unit step), and a = 5 (an exponential impuls), compare the first two columns in Tab. 4. 

The resultant voltage responses computed according to (49) are shown in Fig. 15, including 
relative errors, where also dependences on Volterra series orders are presented. 

The relative errors above were computed via a Matlab ODE45 Runge-Kutta function applied 
directly to the nonlinear ODE (44). As expected, more Volterra terms lead to more accurate 
results, see also (Brancik, 2009) where up to 2"4-order terms were considered, and respective 
Matlab listings are presented. 
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Fig. 14. Numerical inversions leading to voltage response Volterra series terms 
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Fig. 15. Resultant voltage responses and relative errors 


4. Conclusion 


The paper dealt with a specific class of techniques for the numerical inversion of Laplace 
transforms, namely based on a complex Fourier series approximation, and connected with 
a quotient-difference algorithm to accelerate the convergence of infinite series arising in 
the approximate formulae. The 1D to 3D NILT techniques have been programmed in the 
Matlab language (R2007b), and most important ones provided as the Matlab function 
listings. To guide readers all the algorithms were explained on selected examples from 
field of electrical engineering, including right callings of the functions. In contrast to most 
others the NILT methods here developed are utilizable to numerically invert complex 
Laplace transforms, leading to complex originals, which can be useful for some special 
purposes. As has resulted from error analyses the accuracies range relative errors from 
10-8 to 10-10 without difficulties which is acceptable for most of practical needs. Based on 
Matlab functions presented, one could further generalize a vector version of the 1D NILT 
function towards a matrix one, enabling e.g. to simulate multiconductor transmission line 
systems, as is shown in (Brančík, 1999), where, however, an alternative technique, so- 
called £ algorithm, has been applied to accelerate the convergence of infinite series. 
According to the author’s knowledge, the paper presented ranks among few summary 
works describing multidimensional NILT techniques, covering Matlab listings beyond, 
based just on a complex Fourier series approximation, and in conjunction with the 
quotient-difference algorithm, which seems to be more numerically stable compared to 
the € algorithm mentioned above. 
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